Nonlinear r-modes in Rapidly Rotating Relativistic Stars 
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The r-mode instability in rotating relativistic stars has been shown recently to have important 
astrophysical implications (including the emission of detectable gravitational radiation, the expla- 
nation of the initial spins of young neutron stars and the spin-distribution of millisecond pulsars 
and the explanation of one type of gamma-ray bursts), provided that r-modes are not saturated at 
low amplitudes by nonlinear effects or by dissipative mechanisms. Here, we present the first study 
of nonlinear r-modes in isentropic, rapidly rotating relativistic stars, via 3-D general-relativistic hy- 
drodynamical evolutions. Our numerical simulations show that (1) on dynamical timescales, there 
is no strong nonlinear coupling of r-modes to other modes at amplitudes of order one - unless 
nonlinear saturation occurs on longer timescales, the maximum r-mode amplitude is of order unity 
(i.e., the velocity perturbation is of the same order as the rotational velocity at the equator). An 
absolute upper limit on the amplitude (relevant, perhaps, for the most rapidly rotating stars) is set 
by causality. (2) r-modes and inertial modes in isentropic stars are predominantly discrete modes 
and possible associated continuous parts were not identified in our simulations. (3) In addition, the 
kinematical drift associated with r-modes, recently found by Rezzolla, Lamb and Shapiro (2000), 
appears to be present in our simulations, but an unambiguous confirmation requires more precise 
initial data. We discuss the implications of our findings for the detectability of gravitational waves 
from the r-mode instability. 

PACS numbers: 95.30.Sf, 04.40.Dg, 97.10.Kc, 95.30.Lz 



Considerable effort has been spent in the last two years, 
in determining the properties of r-modes |IJ in rotat- 
ing compact stars, since the discovery that these 
modes are unstable to the emission of gravitational ra- 
diation. This is motivated by the current understanding 
that the r-mode instability may have several important 
astrophysical consequences: it provides an explanation 
for the spin-down of rapidly rotating proto-neutron stars 
to Crab-like spin-periods and for the spin-distribution of 
millisecond pulsars and accreting neutron stars, while be- 
ing a strong source of detectable continuous gravitational 
radiation (see for reviews). In addition, if r-modes 
induce differential rotation, then the interaction between 
them and the magnetic field in neutron stars has been 
proposed as a gamma-ray burst model [Q or as a mech- 
anism for enhancing the star's toroidal magnetic field, 
which could, in turn, limit the r-mode amplitude ||. 

Still, several important uncertainties remain, the reso- 
lution of which could significantly modify the above con- 
clusions. In this Letter we present nonlinear hydrody- 
namical simulations of r-modes in rapidly rotating rela- 
tivistic stars (in the relativistic Cowling approximation 
H ) and address the questions of the maximum amplitude 
that r-modes can reach, the nature of their frequency 
spectrum, and the existence of a kinematical differen- 
tial drift, associated with r-mode oscillations ||. To our 
knowledge, these are the first simulations of nonaxisym- 
metric oscillation modes in rotating (Newtonian or rela- 
tivistic) stars. 

The maximum amplitude of unstable r-modes in a fluid 
star (neglecting the magnetic field) and the precise mech- 



anism by which it is set, are currently unknown, but sat- 
uration of the amplitude is thought to occur due to some 
form of nonlinear hydrodynamical coupling. An exam- 
ple of such a mechanism is the nonlinear coupling of the 
unstable r-mode to other, stable, modes of pulsation. 
Presumably, nonlinear saturation is set on a hydrody- 
namical timescale, although it cannot be excluded that 
weak hydrodynamical couplings saturate the r-mode am- 
plitude on longer timescales (but shorter than the growth 
timescale due to gravitational radiation reaction). 

For our study we have used a numerical code 
based on the 3-D CACTUS code developed by the 
AEI-Potsdam/NCS A/ Washington University collabora- 
tion |Io| , in which we implemented the 3rd order Piece- 
wise Parabolic Method (PPM) [jllj for the hydrodynam- 
ics (see [^2| for a recent review) and added initial data 
for equilibrium and perturbed rapidly rotating relativis- 
tic stars. In |ll| it was shown that the 3rd order PPM 
method is suitable for long-term evolutions of rotating 
relativistic stars. The equilibrium initial data are con- 
structed using the RNS code We focus on a par- 
ticular, representative, rapidly (and uniformly) rotating 
model with gravitational mass M = 1.63M Q , equato- 
rial circumferential radius R — 17.25km and spin pe- 
riod P = 1.26ms (the gravitational mass and equatorial 
radius are larger than for a corresponding nonrotating 
model, because of rotational effects - the ratio of polar to 
equatorial coordinate radii is 0.7). The star is rotating at 
92% of the mass-shedding (Kepler) limit (at same central 
density) and we use the N = 1.0 relativistic polytropic 
equation of state (see W). Unless otherwise noted, we 
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use 116 3 Cartesian grid-points, with an equidistant grid 
separation of Ax = Ay = Az = 0.31km. 

During the time-evolution, we only evolve the hydro- 
dynamical variables, keeping all spacetime variables fixed 
at their initial, unperturbed, values (for small ampli- 
tude pulsations this is equivalent to the Cowling approx- 
imation in linear perturbation theory; see ]l3]]). Since 
r-modes are basically fluid modes, we expect that the 
adopted approximation is suitable for studying, quali- 
tatively, their basic properties. The computational re- 
quirements for a coupled spacetime and hydrodynamical 
evolution (for the same long-term accuracy as reported 
here using the Cowling approximation), by far exceed 
presently available supercomputing resources. 
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FIG. 1. Time-evolution of the rotational velocity profile 
for a stationary equilibrium model, using the 3rd order PPM 
scheme. The star remains stationary, even after 20 rotations. 
Larger deviations (due to finite-differencing) occur only near 
the surface of the star. 

Figure 1 clearly demonstrates the ability of the 3-D 
code to maintain an unperturbed stationary equilibrium 
configuration (apart from the influence of local trun- 
cation errors) in long-term evolutions. The rotational 
velocity profile (v x along y-axis, where are the con- 
travariant components of the equilibrium 3-velocity, as 
measured by a local zero-angular-momentum observer) 
remains nearly unchanged over more than 20 rotational 
periods (all simulations presented here are stopped due to 
computing time restrictions only and could be continued 
for longer times). 

In order to excite r-modes during the time-evolution, 
we perturb the initial stationary model, by adding a spe- 
cific perturbation 8v l to the equilibrium 3-velocity v l . 
As there is no exact linear eigenfunction available in the 
literature that would correspond to an I = m = 2 r- 
mode eigenfunction for rapidly rotating relativistic stars 
in the Cowling approximation, we are using an approx- 
imate eigenfunction, derived in spherical polar coordi- 
nates, valid in the slow-rotation Q(£l) limit to the first 
Post-Newtonian (1PN) order (see |ll|[l6| for details). We 
apply the usual transformation between the spherical po- 
lar and Cartesian coordinate systems as the latter is used 



in the 3-D code. We note that, in the Newtonian limit, 
an amplitude of a = 1.0, corresponds to a maximum 
velocity perturbation (i.e. the maximum value of the ve- 
locity component Sv e in the equatorial plane) equal to 
roughly 1/3 of the rotational velocity of the star at the 
equator. With this eigenfunction (multiplied by a cho- 
sen initial amplitude a, that coincides, in the Newtonian 
limit, with the definition of the amplitude in [Jf^]), we 
are able to excite mainly the I — m = 2 r-mode. Due to 
the approximate nature of the eigenfunction, additional 
modes are also excited, primarily m = 2 inertial modes. 

Figure || displays the evolution of the axial velocity 
in the equatorial plane (v z along the y-axis) at a coor- 
dinate radius of r — 0.75r e , where r e is the coordinate 
radius at the equator. It is clear that the evolution is a 
superposition of several modes, and that one mode, the 
I = m = 2 r-mode, is the dominant component. The 
amplitude in this evolution is a = 1.0. The perturbed 
star is evolved for more than 25ms (26 I = m — 2 r-mode 
periods), during which the amplitude of the oscillation 
decreases due to numerical viscosity. Even at amplitudes 
larger than 1.0, the evolution is still similar to that in 
Figure |^, with no sign of nonlinear saturation of the r- 
mode amplitude on a dynamical timescale (only when 
a exceeds a value much larger than 1.0, nonlinear hy- 
drodynamical saturation sets in; however the precise de- 
termination of the maximum saturation amplitude will 
require more accurate inital data and the simultaneous 
dynamical evolution of the gravitational field). Our re- 
sult implies, that, unless nonlinear saturation is set at 
timescales much longer than the dynamical one, nonlin- 
ear hydrodynamical couplings would not prohibit gravi- 
tational radiation reaction to drive unstable r-modes to 
large amplitudes (of order one) before saturation sets in. 
An absolute upper limit on the r-mode amplitude (rele- 
vant, perhaps, for the most rapidly rotating stars) is set 
by causality, requiring \J ViV 1 < c (where c is the speed 
of light), or approximately, a C ausai ^ 3c/QR, where fi is 
the angular velocity of the equilibrium star (for slowly 
rotating stars, we expect other upper limits to be more 
relevant). 

A Fourier transform of the time-evolution shown in 
Figure ||, as a function of the frequency in the inertial 
frame, shown in Figure |^, reveals that the initial data we 
are using, excite mainly the I = m = 2 r-mode (r-z), with 
a frequency of 1.03 kHz and, with smaller amplitudes, 
several inertial modes (13, 14, 15) and higher harmonics. 
The frequencies of the various modes are the same at any 
given point inside the star (the specific data shown in the 
figure correspond to radii r = 0.5r e and r = 0.75r e ), i.e. 
the evolution consists of a sum of discrete modes and pos- 
sible associated continuous parts were not excited. This 
is in agreement with the conclusions in p6| , regarding 
the r-mode spectrum in isentropic stars. The possible 
existence of a continuous part in the frequency spectrum 
of nonisentropic stars (see e.g. pq]), will be examined in 
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FIG. 2. Evolution of the axial velocity in the equatorial 
plane for an amplitude of a = 1.0, at r = 0.75r e . The evolu- 
tion is a superposition of (mainly) the I — m — 2 r-mode and 
several inertial modes. The amplitude of the oscillation de- 
creases due to numerical (finite-differencing) viscosity of the 
code. A beating between the I = m = 2 r-mode and the 
lo — 4, m = 2 inertial mode can also be seen. 

In order to identify some of the peaks in the Fourier 
transform with specific inertial modes, we compare the 
ratio of their frequencies over the frequency of the I = 
m = 2 r-mode to the corresponding ratios derived from 
the normal mode linear eigenfrequencies, for a Newto- 
nian N = 1.0 polytrope JT^j . This comparison shows 
good agreement and allows us to identify, in Figure [|, 
the m = 2 inertial modes with l a = 3,4,5 Jl9), having 
corresponding frequencies / = 0.68, 1.16 and 0.38 kHz. 
In the Newtonian limit, one would expect the above in- 
ertial modes to have frequencies / = 0.70, 1.14 and 0.37 
kHz, for an I = m = 2 r-mode frequency of 1.03 kHz, 
which shows that, even though the effects of relativity 
and rapid rotation on the individual frequencies is not 
negligible, the ratios of the r-mode frequency to the fre- 
quencies of the inertial modes remain close to their New- 
tonian, slow-rotation values. 

We investigate the influence of the numerical viscosity 
of the code (due to finite-differencing) on the evolution, 
by comparing evolutions at different grid-spacings. Fig- 
ure [I] shows the late-time evolution of the axial velocity 
in the equatorial plane, at r = 0.75r e for grid-spacings 
of Ax = 0.59km and Ax = 0.31km. The decrease in the 
amplitude of the oscillation (when compared to the initial 
amplitude) scales as roughly first order with grid-spacing. 
This indicates that the numerical viscosity damping the 
oscillations is dominated by the truncation error at the 
surface of the star, which is only first-order (compared 
to the second-order accuracy in the interior; see related 
discussion in pjfl). 
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FIG. 3. Fourier transform of the evolution shown in Figure 
^[ showing the frequencies of the modes in the inertial frame. 
The frequencies of the modes are the same at different radii 
(discrete spectrum) and the ratio of the I = m = 2 r-mode 
frequency to the frequencies of the identified inertial modes 
agrees with published results. 

In a recent paper || it is indicated that the hydro- 
dynamical oscillatory motion of an r-mode in a rotating 
star may be accompanied by a differential drift (of kinc- 
matical origin). In the equatorial plane, the drift velocity 
has opposite sign, compared to the rotational velocity of 
the unperturbed star. In a rotating star with a poloidal 
magnetic field, this kinematical drift may wind up the 
magnetic field lines and limit the effect of the r-mode 
instability The differential drift reported in || is 

second order in the r-mode amplitude, but it is derived 
from the first-order in amplitude definition of the r-mode 
velocity. Still, in the restricted case of a spherical shell 
p0| , the full nonlinear hydrodynamical equations yield 
the same kinematical drift as in || , independent of the 
amplitude of the oscillations. 

Using our code we attempted to investigate the pres- 
ence of such a kinematical drift in the numerical evolu- 
tions. Our numerical results do indicate that the per- 
turbed star is rotating slower near the surface (compared 
to the rotational velocity of the unperturbed star) and 
that this drift scales roughly as a 2 for amplitudes of order 
a ~ 1.0, although its magnitude is significantly smaller 
than estimated in ||. We cannot test the presence of 
the drift for amplitudes much less than a ~ 1.0, since it 
becomes smaller than numerical truncation errors. Even 
at a ~ 1.0, it is not completely unambiguous that the 
drift we find is the same as that predicted by as it 
could be due to the fact that our initial data do not 
correspond exactly to a single I = m = 2 r-mode and 
thus also contain other modes and violate (to some ex- 
tent, which is still acceptable for the main purpose of this 
work) the relativistic Hamiltonian and momentum con- 
straints. However, the fact that our drift scales as a 2 , 
points to an association with the results reported in [||. 
Our preliminary findings on the differential drift need to 
be confirmed by simulations with more precise eigenfunc- 
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tions and higher resolution, before a definite conclusion 
can be drawn. 
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FIG. 4. Late-time evolution of the axial velocity in the 
equatorial plane for a = 1.0, at r = 0.75r e , for two differ- 
ent grid-spacings. The decrease in the oscillation amplitude 
scales as roughly first order with resolution (it is dominated 
by the first-order truncation error at the surface of the star). 

In summary, our hydrodynamical evolutions of non- 
linear r- modes show that (1) nonlinear hydrodynami- 
cal couplings would not prohibit r-modcs in isentropic, 
rapidly rotating relativistic stars from attaining a large 
amplitude, of order one, when driven unstable by gravi- 
tational radiation reaction - it remains to be investigated 
whether nonlinear saturation could set in at timescales 
longer than the dynamical one, (2) r-modes and inertial 
modes, in isentropic stars, are discrete and (3) a kinemat- 
ical drift associated with r-mode oscillations, appears to 
be present in our simulations, although an unambiguous 
confirmation will require more precise initial data. 

Our finding that gravitational radiation could drive 
unstable r-modes to a large amplitude, implies that r- 
modes can easily melt the crust in newly-born neutron 
stars leaving the initial conclusions about the r- 

mode instability being a strong source of gravitational 
waves (see ||||) essentially unchanged. Our results also 
imply that r-modes could be excited to large amplitudes 
during the violent formation of a proto-neutron star, after 
a supernova core-collapse ^] . Our finding that r-modes 
oscillate at a, predominatly, discrete frequency (at least 
in isentropic stars) significantly simplifies the expected 
gravitational wave signal compared to a signal emitted 
by a source with a significant continuous part in its fre- 
quency spectrum. The indications for the existence of 
the kinematical drift, call for a more detailed consider- 
ation of the interaction between r-modcs and magnetic 
fields. 

We will present details and tests of our numerical 
method, as well as extensive higher-resolution results for 
various equations of state and rotation rates, in a forth- 
coming paper ]23| ]. In future work, we plan to study 
r-modes in nonisentropic stars and implement an accel- 



erated gravitational-radiation reaction force (see p4| ) . 
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